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Abstract 

This paper is a companion to a GDC 2008 Lecture with the same title. It provides a brief 
overview of spherical harmonics (SH) and discusses several ways they can be used in interactive 
graphics and problems that might arise. In particular it focuses on the following issues: How to 
evaluate lighting models efficiently using SH, what "ringing" is and what you can do about it, 
efficient evaluation of SH products and where they might be used. The most up to date version 
is available on the web at http://www.ppsloan.org/publications 


Introduction 

Harmonic functions [2], the solutions to Laplace's equation, are used extensively in various 
fields. Spherical Harmonics are the solutions when restricted to the sphere’. They have been 
used to solve potential problems in physics, such as the heat equation (modeling the variation of 
temperature over time [5][25]), and the gravitational and electric fields[9]. They have also been 
used in quantum chemistry and physics to model the electron configuration in atoms and model 
quantum angular momentum [16][51]. Closer to graphics they have been used to model 
scattering phenomena [7][17]. In computer graphics they have been extensively used, early 
uses were in modeling volumetric scattering effects [18], environmental reflections for micro- 
facet BRDF's without global shadows[6], non-diffuse off-line light transport simulations[40], 
BRDF representations [53], image relighting[28], image based rendering with controllable 
lighting [54][55], and modeling light source emission[8]. More recent examples include more 
work in atmospheric scattering [50] and computer vision [3]. 


The focus of this article is on techniques related to interactive rendering. The first paper that 
has been used extensively in games deals with using Spherical Harmonics to represent 
irradiance environment maps efficiently, allowing for interactive rendering of diffuse objects 
under distant illumination [35]. This was extended to handle a limited class of BRDF's with the 
same constraints [36]. Precomputed radiance transfer (PRT) [41][20][24] models a static 
object/scenes response to a lighting environment, often represented using SH, including 
complex global illumination effects like soft-shadows and inter-reflections with diffuse and 
simple glossy materials. This was extended to handle more general BRDF models [20][23] [42], 


Т They are the angular portion of the solution in spherical coordinates. 


incorporate sub-surface scattering [42], significantly improve rendering efficiency through 
compression techniques from machine learning [42], and various techniques to model “local” 
texture like surface details [43][44][45]. SH have been used to model single scattering from 
distant lighting environments [49]. Other uses have been in using gradients to lift the distant 
lighting assumption [1], several techniques addressing dynamic objects [56][37] including 
support for inter-reflections [46][33], as a representation of visibility to model shadowing of 
objects with general BRDF models [12], using scaling operators to model shadows from 
deformable objects [52], as a parameterization of refraction [11] and as a technique to address 
the level of detail problem with normal maps [15]. 


More practical papers include covering implementation details for PRT [13], how to integrate 
these techniques into an engine [30], how to use SH+gradients for irradiance volumes [31][32], 
practical issues around projection and how to efficiently quantize SH coefficients [21] and a nice 
paper that projects an analytic skylight model [34] into SH and uses a global polynomial fit to 
evaluate the SH lightprobe as a function of the models parameters [14]. Numerical techniques 
for more robustly projection functions defined over the hemisphere using SH have also been 
investigated [22]. 


Many of the uses in real-time graphics are just as a convenient way for representing spherical 
functions — visibility, lighting and reflectance. While there are many other basis functions that 
can be used, wavelets [39], wavelets on cube maps [27], spherical radial basis functions [9], and 
others [26], spherical harmonics have some nice properties that will be described in this 
document. It is important to stress that there are scenarios where these other basis functions 
are more appropriate. 


While spherical harmonics may seem somewhat daunting, they are actually straightforward. 
They are the spherical analog to the Fourier basis on the unit circle, and are easy to evaluate 
numerically. Like the Fourier basis used in signal processing, care has to be taken when 
truncating the series (which will always be done in video games), to minimize the "ringing" 
artifacts that can occur. This article will describe how to evaluate and represent lights efficiently 
using spherical harmonics, how to pull conventional lights out of a SH representation, describe 
"ringing" and mitigation techniques to minimize its impact, and go over products of functions 
using spherical harmonics, describing where they are useful and special cases that are worth 
optimizing for. 


Background 
Definition Spherical Harmonics define an orthonormal? basis over the sphere, S. Using the 
parameterization 


5 = (x,y,z) = (sin соѕф,ѕіп Ө sing, cos 0) 


2 An orthonormal basis is one that has the property S fif; = бу. that is one if i==/, and zero otherwise. 


Where s are simply locations on the unit sphere. The basis functions are defined as 
Y" (8,9) = КГ ene B" (cos 6), LEN,-l<m<l 


Where Р)" are the associated Legendre polynomials and Kj” are the normalization constants 


(21+ 1)(0 — т]! 
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The above definition is for the complex form (most commonly used in the non-graphics 
literature), a real valued basis is given by the transformation 


V2Re(Y/*) m > 0 V2K/" cosmo Pi" (cos8) т > 0 
yi” = AJ2Im(Y7) m < 0 = 4/2kKP" sin[m|g Р," (cos 0)т «0 
үд т-0 KPP? (cos 9) m=0 


The index / represents the “band”. Each band is equivalent to polynomials of that degree (so 
zero is just a constant function, 1 is linear, etc.) and there are 2/+1 functions in a given band. 
While spherical coordinates are convenient when computing integrals, they can also be 
represented using polynomials, as is commonly done when evaluating them (see Appendix A1 
Recursive Rules for Evaluating SH Basis Functions and Appendix A2 Polynomial Forms of SH Basis 
for details.) An order n SH expansion uses all of the basis functions through degree n-1. 


Spherical Harmonics can be visualized т a couple of ways. One standard way is to distort a unit 
sphere, by scaling each point radialy by the absolute value of the function and coloring it based 
on the sign (red for positive, blue for negative.) Above are images of the first three bands using 
this technique. 


The functions in the central column (|=0) are called zonal harmonics (ZH) and will be discussed 
later, these functions have rotational symmetry around the z axis and the zeros (locations where 
the function is zero) are contours on the sphere parallel to the XY plane. The functions where 
(I ||) are called sectorial harmonics and the zeros define regions like apple slices. 


An alternative visualization is to draw them using the parameterization of a cube map unfolded 
onto the plane. The unfolding of the cube map is as follows: 


Here magnitude is encoded with color (red positive, blue negative, zero green) and iso-intensity 
contours have been evenly placed (white lines) to give more intuition for the gradient of the 


function (when they bunch together the function is changing faster, etc.) 


T 
T 


Projection and Reconstruction Because the SH basis is orthonormal the least squares 
projection of a scalar function f defined over S is done by simply integrating the function you 
want to project, f (s), against the basis functions (proof in Appendix Аб Least Squares 
Projection) 


fn = | foro ds 


These coefficients can be used to reconstruct an approximation of the function f 


n l 
f=) > f" 
1-0 т--1 
Which is increasingly accurate as the number of bands n increases. This paper concentrates оп 
low-frequency approximations to f, for higher frequency representations other bases tend to do 


a better job. Projection to n-th order generates n° coefficients. It is often convenient to use а 
single index for both the projection coefficients and the basis function, via 
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f(s) = У fils) 
1-0 


Where i=/(/+1)+m . This formulation makes it clear that evaluating at direction s of the 
approximate function is simply a dot product between the n° coefficient vector f; and the vector 
of evaluated basis functions y;(s). The first coefficient (fy or fo depending on indexing) 
represents the average value of the function over the sphere and will sometimes be referred to 
as the DC term. 


Basic Properties One important property of SH is how projections interact with rotations. 

Given a function g(s), which represents a function f(s) rotated by a rotation matrix Q, so g(s) = 
f(Q(s)) the projection of g is identical to rotating f and re-projecting it. This rotational 
invariance is similar to the translational invariance in the Fourier transform. This means that, for 
example, lighting will be stable under rotations, so there won't be any aliasing artifacts or 
“wobbling” of the light sources. Below are images of a sphere illuminated by a directional light 
source, the top row is using SH, the bottom row is using the Ambient Cube basis? from Valve 
[26]. The first column is a best case orientation, and the second is near a worst case one. The 
image is invariant using SH. This basis is discussed in more details in Appendix A9 Ambient Cube 
Basis. This will happen to some extent with any other basis defined over the sphere. 


? This basis is more efficient to evaluate compared to SH, having 6 basis functions with only 3 being non- 
zero at any point on the sphere. 
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Due to the orthonormality of the SH basis, given any two SH functions a and b, the integral of 


the product is simply the dot product of the coefficient vectors: f ü(s)b(s) = " aibi 


Convolution Given a kernel function h(z) that has circular symmetry, you can generate a new SH 
function that is the result of convolving the kernel with an original function f. h must have 
circular symmetry for the result of the convolution to also be represented on the sphere S, 
instead of the rotation group SO(3). The convolution can be done directly in the frequency 
domain using the following equation: 
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(h *f)r = air 


This amounts to simply scaling each band of f by the corresponding m=0 term from h. 


Rotation As mentioned before, the SH are closed under rotation. SH rotation matrices are in a 
block structure, where each band is rotationally independent and has a dense (21+1)х(21+1) sub- 
matrix. There are several ways to compute these rotation matrices, for very small orders 
(quadratic and less) doing so symbolically is most efficient, but for higher orders it seems to be 
more efficient to decompose the rotation matrix into zyz Euler angles [19]. 


Zonal Harmonics Spherical harmonic projections of functions that have rotational symmetry 
around an axis are called Zonal Harmonics (ZH.) If they are oriented so that this axis is Z, the 


zeros of the function will form lines of constant latitude and the functions only depend on 0. 
The coefficient vector in this orientation only has one non-zero per band, so a n-th order 
function has n instead of n^ coefficients. Zonal Harmonics have been used to approximate 
transfer [44], and are the common representation of phase functions in scattering theory 
171117 |, they will be used extensively in this paper when modeling light sources. Rotation of 
Zonal Harmonics is simpler than general SH, it can be done with what is effectively a diagonal 
matrix and only requires evaluating the SH basis functions in the new direction d. Given the ZH 
coefficients of a function (only the т=0 terms from an SH projection) z; it can be rotated to a 
new direction d using this equation: 


=> a ale (ауу) 
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So the resulting SH coefficients are: 
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4r 
21-1 
SH Products The К" coefficient of the product of two functions f and g represented using n^ 
order SH projected into SH has the following form: 


pe = | е Улиг palay nets 
i=0 150 Ü 


WhereT is the triple product tensor: 


гу = [rOy Oas 


an order 3 sparse symmetric tensor. Since SH are polynomials, the polynomial product уу have 
maximal degree 2n-2, which means it will have non-zero coefficients through order 2n-1. This 
becomes unwieldy as the number of functions being multiplied grows, so it is common to 
truncate the product early [56][37]. The number of non-zero coefficients as a function of n is 
quite large [47][37] so care has to be taken when generating efficient code. One special case 
that is useful to point out is that if the function f is fixed (ie: distant lighting) you can compute a 
"product matrix" M; which will significantly reduce the cost. This matrix is symmetric and built 
using the following equation: 


(Mp) = 2 Tijk fk 
k 


Computing the product with a function g in this case is simply a matrix vector product. 


Irradiance Environment Maps 

An irradiance environment map is created by convolving a light probe with a clamped cosine 
function; this should be normalized by dividing by z to display radiance. This convolution can be 
done efficiently using SH [35], and is accurate enough to be efficiently rendered directly from SH 
as well. Order 3 SH do a good job approximating this kernel, but if HDR light sources are going 
to be used you might want to consider using Order 5 (the order 4 ZH coefficient is zero so that 
band can be skipped.) Below are images of the clamped cosine kernel and the order 3 SH 
approximation, the red curve is the SH approximation, the figure on the left is a plot as a 
function of theta, on the right a polar plot scaled by the absolute value of the function: 


The Order 3 SH approximation over estimates by 1/16" at theta-0 (north pole) and has a 
spurious lobe at the south pole with magnitude of 1/16". A directional light source that would 
reflect a value of 17 when a normal points at it would reflect a value of 1 pointing in the 
opposite direction (should be reflecting O.) The order 5 approximation has a negative lobe, 
which would reflect -1 with a directional light that would reflect 31 with a normal pointing right 
atit. While these approximations are accurate, the approximation can cause error, particularly 
if very bright light sources are being used. 


Appendix A10 Shader/CPU code for Irradiance Environment Maps contains shader and CPU code 
for efficiently evaluating irradiance environment maps. 


Lighting Models 

There are various ways of representing lighting in SH. The simplest is to just project from a cube 
map, but there are also analytic models that are inexpensive to evaluate and potentially useful 
to expose to artists. A recent paper [14] does a nice job of projecting a practical skylight model 
[34] into SH and fits a global polynomial of the SH coefficients over the parameter space of the 
model. 


Projection from Cube Maps 

To project from a cube map you simply need to integrate the SH basis functions against the cube 
map. This can be done numerically by evaluate the SH basis functions in the direction of each 
texel center, weight it by the differential solid angle for that texel and normalize the results. In 
pseudo-code: 


float f[],s[]; 
float fwtSum=0; 
Foreach(cube map face) 
Foreach(texel) 
float fTmp = 1 + u^2-*4v^2; 
float fwt = 4/(sqrt(fTmp)*fTmp); 
EvalSHBasis(texel,s); 
f += t(texel)*fWt*s; // vector 
fwtsum += fWt; 
f *= 4*Pi/fwtSum; // area of sphere 


In the code above и and v represent the two coordinates on the given face that are not +/- 1°, 
t(texel) is the texel color. EvalSHBasis needs to normalize the input coordinate (on a cube face), 
and simply evaluate the SH basis functions in that direction. The last normalization can be 
omitted (instead you can just divide by the number of samples), since the normalized sum of the 


^ le: on the +X face, it would be y/z, etc. 


differential solid angles should be 4*Pi, but it tends to be a little ОН (particularly when using low 
resolution cube maps.) 


Below are images of the reconstruction of a HDR light probe into order 1 to 6 Spherical 


Harmonics. The final image is the light probe that was projected. 


Analytic Models 
Directional lights are trivial to compute, you simply evaluate the SH basis functions in the given 
direction and scale appropriately (see Normalization section.) Spherical Light sources can be 
efficiently evaluated using zonal harmonics. Below is a diagram showing an example scene, we 
want to compute the incident radiance, a spherical function, 
at the receiver point P. Given a spherical light source with 
center C, radius r, what is the radiance arriving at a point P d 
units away? The sin of the half-angle subtended by the light 
source is r/d, so you just need to compute a light source that 
subtends an appropriate part of the sphere. The ZH 
d coefficients can be computed in closed form as a function of 


this angle: 2; = 1: 1 уд dod6 where a is the half-angle 


subtended. See Appendix A3 ZH Coefficients for Spherical 


Light Source for the expressions through order 6. 
Figure 1 Spherical Light Source 


The spherical technique can also be used to model а cone (think of it as a disc at infinity) with 
constant emission. A softer cone can be modeled that has a smooth fall-off over the visible 
portion? — see Appendix A4 ZH Coefficients for Smooth Cone for the equations. 


The columns are for orders 4 and 6, the top row is the projection of a cone and the bottom row 
is a cone that has a smooth fall off. The angles are 90 (green), 45 (red), 30 (blue), 12.5 (black) 
degrees. The dashed lines represent the actual function (slightly displaced for the cone so they 
don't overlap), the solid lines are the SH approximation. In general the soft cone is better 
behaved. How to address the artifacts that arise in projection is the subject of the section on 


Ringing below. 
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Order 6 


? An "ease spline" is used, this is a cubic polynomial f(x) that has the constraints 
f(0)-1,f (0)=0,f(a)=0,f (a)=0, the unique cubic polynomial is: f (x) = 2%° в - ЗА +1 


Normalization 

If [0,1] lighting is being used, it is convenient to normalize the radiance vector so that the 
reflected radiance for an unshadowed receiver with a normal pointing directly at the light 
source would be 1.0. Mathematically you want to compute a scale factor c, that when 
multiplied by the lighting vector L will result in unit reflected radiance when integrated against a 
vector T that represents an unoccluded clamped cosine lobe (projection of normalized clamped 
cosine into SH.) So you have: 


1 | 18) T(s)d yields 1 
= —— ГЕ РНН 

с1 (5) T(s)ds C dot(L,V) 
Only the bands that are going to be used for rendering should be used when computing the 
normalization factor. Aligning T with +Z gives a simple analytic formula; here are the 
coefficients of T for the first 6 bands: 


1 v3 45 А —1 
2/n З/п’ NVT ’16/л’ 


For analytic lights you can use an analytic normalization term, Гог the cone light of angle a that 
would be: 


1 
sin? а 
However the resulting lights will not reflect unit radiance, because the projection error of both 
the clamped cosine function and the light source will not be taken into account. For directional 


А m . 16 i ЭТР кет 
lights? the normalization factor is > and for “ambient” lights it is 27r. 


Extracting Conventional Lights from SH 

Given an SH lighting vector, it is possible to approximate it as a single directional light source 
and an ambient light source. This has been used on hardware that does not support vertex 
shaders. Mathematically we want to compute the intensity of a directional light (c) and the 
intensity of an ambient light source (a), so that squared error of reflected radiance is minimized 
for any surface normal (N). Assuming a fixed direction (d) for the light source, the error function 
we want to minimize is: 
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E(c, a) = | (emo +a- f n comic) ds) ам 


Where Hy(4) = max(dot(N, х) /т, 0) is the normalized clamped cosine. If the lighting Le (5) is 
expressed in SH, this has a simple solution. The irradiance environment map represented by the 
new lighting environment should be as close to the input irradiance environment map — 


ҒАН : : "T 2 . 32 
5 This is assuming no higher than 4" order lighting, for 5" or 6" order the normalization factor is == 


minimizing squared error of these two environment maps is equivalent to minimizing reflected 
radiance over all normals: 


E(e,a) = | (cha(s) ниб) + aa G) Нуб) — Le(s) * Hy G2) ds 


Here La(s) is a normalized SH directional light in direction d, and L, (s) is a normalized SH 
constant light (just depends on DC.) The optimal values a and c are: 


7 ЖЕН 
с = 3167 doto (Ly, Ге) 


The lighting vectors above are all turned into irradiance environment maps through convolving 
with the normalized clamped cosine kernel. The dot product above ignores the DC term and 

Ж [0] is the DC term of the lighting environment. The above assumes а known direction, а good 
candidate direction is the "optimal linear" direction [44], formed by normalizing the vector 
(=L. [3], ^£, [1], e [2]), which are the linear coefficients of the lighting environment. 


Extracting Multiple Lights 

It is also possible to extract multiple lights from a SH light probe. This might be done to combat 
ringing (analytic lights won't have negative lobes), to model glossy reflection (just from the light 
sources pulled out in this manner) or to use a small number of shadow zbuffers (pull the lights 
out for both diffuse and glossy parts of the BRDF.) The optimal linear direction works well when 
a light probe is dominated by a single light source, when this isn't the case the directions and 
intensities of the lights should be optimized for in some fashion. One thing to do is to climb "up 
hill" to find local maxima of the function. Given how smooth SH are, for a given order there is a 
finite distance where any point that distance from a distinct peak is guaranteed to reach the 
point using gradient ascent. A set of points can be generated with the property that the point 
on the sphere farthest from the center of its Voronoi’ cell is less then this distance. If you start 
your search from each of these points you should find all of the point local maxima. These 
distances can be found by looking at the projection of a delta function and computing the 
angular distance between peaks and the zeros. Using a conservative estimate of two thirds this 
radius the number of points needed per order are {1,3,6,10,15,22} for the first 6 orders. 


This point set can be computed by computing a simulation: All points have a force with a certain 
fall-off (say 1/d? where d is the Euclidean distance between them) that acts on every other 
point. Generate the net force acting on each point, and subtract off the normal component (so 
it is just a tangential force), move the points by this force with a weight w. If the net sum of 


7 Given a set of points X on the sphere, the Voronoi cell of a given point x; consists of all the points on the 
sphere that are closer to x; than any other point in the set. 


forces increases, decrease w by half and try again, if it decreases double w and try again. This is 
solving for a set of electrons that minimize the sum of the electrostatic charges.. 


Given a SH lighting environment L(s), you want to find a point on the sphere that is a local 
maximum, which is the same as finding a local minima of —L(s). This is a non-linear optimization 
problem, a small number of BFGS iterations [29] converge on the peaks in my experience. The 
gradients of the basis functions need to be computed when using optimization methods. For a 
point on the sphere differentiating the polynomials is trivial, but you need to allow the point to 
go off the sphere when doing a line search, so you want to use symbolic inputs that normalize 
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the coefficients (ie: қ ;) in place of f (x, y, z)). This is done by 


Мх2+у2+22' (х xy? z2' хачу? 
simply computing the gradient at the normalized location, and multiplying it by the Jacobian of 
the normalization function: 


І2-х2 -ху —х2 
[3 L L 
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-х2 -у2 12-22 
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Where L is /x? + y? + 22. 
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Original 3 direction 2 directions 1 directions 


These images compare radiance (top row) and irradiance (bottom row) lighting environments 
when approximating with 3, 2 or 1 directional lights and a ambient light. Given the N most 
significant peaks (in terms of magnitude) this solves for directional and ambient intensities in a 
least squares sense. 


The equations to compute the light intensities for two lobes are: 


A -В 
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Where: 
А = Lao ° Lao, B = Lao o Lat 


The ambient term is then: 


a = (L,[0] — CcoLao[0] + аши0руут 


For three lights the intensities are: 


A?^—D? CD-—AB BD-AC 
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Where: 
C = Lao ° Laz, D = Lay ° Laz, Е = 2BCD  A(4? — B?-C? — D?) 


This matrix is symmetric, and the ambient coefficient is: 


a = (L,[0] — (coLao[0] + ciLa3[0] + слао) 


The derivations of the equations are in Appendix A5 Solving for Coefficients to Approximate SH 
Environment Map with Directional and Ambient Lights. 


Adding all of the degrees of freedom (directions and intensities) and using a non-linear solver 
would generate higher quality results, possibly using the technique above as an initial guess. 


Ringing 

Ringing, also called Gibbs Phenomenon, is a common problem in signal processing. When a 
signal with a discontinuity is projected into a finite Fourier basis (which can only represent 
continuous functions) overshoot and undershoot will happen around the discontinuity. 
Functions that don't have discontinuities can exhibit similar behavior if the projection is 
truncated. We have already seen these problems when looking at lighting models, and in 
representing irradiance environment maps (projection of clamped cosine function.) A similar 
problem occurs in surface design, where when trying to satisfy a set of geometric constraints 
unwanted oscillations can occur. There are two general solutions to these problems: 


1) Windowing the truncated projection coefficients using sigma factors. This is the most 
common solution in signal processing, and can be trivially used with Spherical 
Harmonics [8][41][37]. 

2) Minimizing some form of variational function (minimizing a measure of curvature for 
example), instead of just the standard least-squares error. This is commonly done in 
computer aided geometric design, but also can be efficiently done using Spherical 
Harmonics [38]. 


Windowing® 

One way to minimize the ringing artifacts is to multiply in the frequency domain (which is a 
convolution in the spatial domain) by a kernel with projection coefficients that taper to zero as 
you approach the cut-off frequency. If this function is a sinc? that is stretched out so that it 
reaches zero at the truncated frequency band, it is called using Lanczos sigma factors”. 
Intuitively what this is doing (in 1D) is convolving in the spatial domain with a tight box function, 
making the function smooth enough to be 
represented without excessive ringing. There are 
more sophisticated ways to attack Gibbs phenomena 
[10][4] but they use the SH coefficients to generate 
piecewise analytic functions that would not be as 
convenient for games. 


In our experience the choice of windowing function is 
not as important as having flexibility to trade off 
between ringing and blurring. The image to the left 


shows the two windowing functions (red is sinc, blue 


is raised cosine lobe — called a Hanning window) 


8 The more common use of the term “windowing” in signal processing is used in the spatial domain — 
when filtering images the spatial extent of the filter is ^windowed", and when taking a FFT of a 1D signal 
the signal might be scaled to make it periodic. In this paper it is being done in the frequency domain. 

? sinc(x) — == where the limit is 1.0 when х-0. 

10 Sometimes the sigma factors will be raised to a power to more aggressively reduce the Gibbs 
phenomena — which is the same as repeatedly convolving the signal. 


scaled for 6" order SH (so they reach zero for the 7" order band, the last value used would be 
evaluated at 5 and the functions would only be evaluated at the integers.) The Hanning function 
decays faster than the Lancosz function, which makes it blur more aggressively. 


Equation W=6 W=10 


(1 + cos =) 


И The Lancosz squared is close to the Hanning function, but decays a bit faster. For the reconstruction of 
the delta function to be a height of 2.25, the Hanning function needs a window of size 12.0105 and the 
Lancosz needs one of size 9.8725. The results look visually indistinguishable. 


Above we show the results using Lanczos and Hanning sigma factors. The signal being projected 
is an order 6 delta function, which is the "peakiest" signal you can project into SH and exhibits 
ringing artifacts. The projection of a delta function is a ZH, so we are showing a cross section of 
the sphere, phi is fixed. Radial magnitude is plotted, the sign of the lobes alternate. 


Looking at all the graphs together (red is the raw projection of the delta function) you can see 
how windowing blurs the signal while eliminating the rings (visible near the origin in the figures.) 


Hanning Lancosz 


Minimizing a Functional 

An alternative approach is to try and minimize some function besides squared approximation 
error. One way to do this is to satisfy a set of constraints (for example exact reconstruction at a 
small number of points) and then use the "slack" variables left over (assuming there are enough 
degrees of freedom) to minimize some error functional[38]. Given the low SH order often used 
in games/graphics, this approach does not seem that practical so | will not spend any more time 
on it. An alternative is to attempt to minimize a norm that penalizes large oscillations. This can 
be done with Spherical Harmonics in a straightforward fashion. The Laplace operator, or 
Laplacian, is the divergence of the gradient of a scalar function; equivalently it is the sum of 
unmixed partial derivatives 


лау? ар? df? 
P= tt ay? а? 


In spherical coordinates on the unit sphere this is 
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The integral of the squared Laplacian is a curvature measure used on the sphere [38]. The 
function we are going to minimize is: 


E(c) = | (f - r6) as +a (мо) as 


We will tweak this a little bit, we already know the raw projection coefficients f", so we want 
to find new coefficients cj" that are as close to the least squares result as possible, while also 
minimizing the weighted squared Laplacian. This can be done in closed form [see Appendix A6 
Least Squares Projection], which given a Л results in the following coefficients. 


c= EM ETE 
l 0A 1)2) 


Note that this amounts to a windowing function that is dependent оп A4. When Л is zero, you 
get the least squares coefficients, when A is infinity you get just the DC term, which has zero 
curvature. One approach to choosing 4 is to solve that it reduces the squared Lalpacian by a 
fixed amount, say one half. This can be done using any standard root finding technique. See 
[Appendix A7 Solving for Lamba to Reduce the Squared Laplacian] for an explanation of how to 
do this using Newtons method. Below are the results with a 6ht order delta function solving for 
the squared Laplacian to be 10% (A = .004209 green) and 50% (A = .000632 blue) of the 
original. The final plot has them along with the delta function itself: 


42=0.1 42-0.5 Composite 


Here are images using an actual lighting environment. The first column of images are contour 
plots (blue is negative) of the second column of images - windowing with a broad (smoother) 
window. The forth column of images are contour plots of the third column — smoothing with a 
narrower (less blurry) filter. The top row is the original image — all results are shown at 6" 
order: 


Original 


Cos 


Sinc 


Windowing should be used judiciously. When using irradiance environment maps, the 
convolution with the clamped cosine function aggressively attenuates high frequencies and 


windowing is rarely necessary. Also scenes where there is a lot of normal variation tend to not 
show off the ringing artifacts as much as ones where the normals of the receiver are smoother. 
Below are images of a simple scene (“доог” and a ground plane) that show how windowing 
affects the shaded result. An area where ringing occurs is highlighted in the second row. 


” 


Мо Windowing Hanning 6 Hanning 10 


When using HDR, or just moderately bright lights, ringing is a more severe problem and even 
irradiance environment maps might need some amount of windowing 


Another problem that can happen with ringing is color artifacts. In the images below we have a 
sphere illuminated with a bright directional yellow light (almost rim lit from the upper right, 
order 6) and a moderate ambient white light. The top row uses irradiance environment maps, 
the bottom row uses order 6 PRT (so a more accurate approximation of the cosine kernel is 
used.) The first column is without windowing, the second column uses a Hanning window of 
order 4 on the top row and order 5 on the bottom. The unwindowed version shows both the 
positive lobe (order 3, negative order 6) and the blue bands are from the negative rings of the 
directional light (remove red/green from the ambient light, but leave the blue alone.) 


Order 3 


Order 6 


No Windowing Hanning 


Content Sensitive Windowing 

Lighting can be globally windowed, but you can also window depending on the influence the 
rings will have on the final shaded image. An example is illustrated below. We are rendering a 
somewhat matte (but not diffuse — phong with power 10) ball using a bright directional light 
source (the highlight is saturating.) If no windowing is applied, the principle highlight is sharper, 
but artifacts from the rings are clear. If the lighting is windowed, the rings disappear but the 
principle highlight is blurred as well. Instead you can look at the angle between the reflection 
vector and the dominant light direction, if this is small, no windowing is necessary, as it 
increases you can blend between the windowed and non-windowed light source, preserving the 
sharp highlight while eliminating the ringing artifacts. This sequence of images is shown below: 


6" order reflection Windowed (Hanning 7) Content sensitive windowing 


The equation used to control the blending is: 


0 (ег) = 2 | 


Wa = тах (о er 
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Where w, is the weight of the un-windowed light source (windowed is 1 — w4), c; is a threshold 
that determines when you completely use the windowed light source and p controls the 
transition region of the blend. In this figure c, is 0.07 and p is 0.8. You can experiment with 
these parameters, the thresholds will largely be based on the amount of windowing and the 
material properties. 


If the light source is not a directional one, you can compute how well it can be approximated by 
a directional light source and take that into account when deciding how much to blend (the 
simplest way would be to compute a tensor product between how close the lighting 
environment is approximated by a directional light and how close the normal/reflection vector is 
to the dominant light vector.) For more complex shading, say for example PRT, the dominant 
transfer direction can be used, and how well it approximates the transfer function can be yet 
another factor in the tensor product (so you window when any one of the terms isn't 
approximated well.) 


If the lighting is very dynamic, you will need to watch out for temporal artifacts that could occur 
as the lights change (ie: shadows/reflections would become sharper, etc.) For things like static 
SH light probes this technique should work well though. 


SH Products 


Computing the SH representation of the product of two functions represented using SH is often 
useful. Example scenarios are: 


3) Punch a hole in a skylight model based on a large flying object (visibility times light), or a 
simple visibility model of the scene (large buildings, etc.) 


4) Multiply visibility functions. This happens when doing dynamic approximate global 
illumination. 

5) Scale or modify a SH light probe. Multiplying by a some constant between zero and 1 
can be used to approximate clouds for example. 


Computing a product in the frequency domain is quite complex, it boils down to applying a 
"triple product tensor" times the two SH vectors. Code for this can be generated efficiently [47] 
and will not be described in this paper. There are several special cases that are worth 
mentioning though. 


Products with a Constant Function 

If one of the SH functions is going to be used a lot you can build a dense matrix called the 
product matrix, this makes the triple product a simple matrix times vector product, which is 
significantly less expensive. An order 6 product would have 1296 multiplies instead of the 2527 
in the code generated by [47]. 


Products with varying Orders 

This is particularly common when the output order is lower, for example quadratic, so that you 
can represent a local radiance environment. Special casing the code for these cases can 
significantly decrease the complexity of the code. For example the product of two 6" order SH 
has 2527/1995 multiplies/adds when computing 6" order results, but only 933/676? when 
computing 3" order results. Another example is simple ambient occlusion, in this case one of 
the terms is just DC, and so you simply have to scale the other vector by the DC value. Finally, 
one of the two functions might be lower order (ie: just multiply by linear visibility) which can 
also reduce the cost. 


Product with Zonal Harmonics 

If one of the functions is a Zonal Harmonic, you can rotate the other function into that same 
frame (less expensive rotation, only need 2 Euler angles due to the symmetry), compute the 
product and rotate back. The sparsity in the ZH frame eliminates a significant amount of work, 
which can increase performance. The product of two 6" order functions, one which is a ZH 
oriented in z, requires only 380/249 multiplies/adds. For arbitrary ZH the time to compute one 
million products (where one of the functions is always a in an arbitrary orientation) is around 1.2 
seconds while the time to compute general 6" order products is just over 3 seconds, so this 
technique is almost 3 times faster. 


Product with Analytic Function 

If one of the functions has an analytic form, it is more accurate to analytically compute what 
amounts to the product matrix. An example is zeroing everything under the horizon (useful if 
there is a ground plane) or taking a product with a clamped cosine function. Doing this 


12 More efficient code could be generated, these results are using a naive algorithm applied to the output 
of [47]. 


analytically is equivalent to having an infinite order expansion of the analytic function, and in 
general will be much faster than using a SH expansion in these cases. The code for these two 
examples is in an appendix [Appendix A8 Code for Multiplying SH by Analytic Functions]. 


Conclusions 

Spherical Harmonics are extremely useful tool, particularly for lighting, in games. Hopefully this 
article shed some light on how they can be used and how to mitigate some challenges that can 
arrive when using them. The are several ways to extend the ideas discussed in this paper, 
windowing coefficients (using any window function) can be solved for that minimize the 
magnitude of a negative (or positive) lobe, or possibly the magnitude of reflected radiance when 
pointing in a direction that should be zero. When extracting lights, more rigorous techniques 
based on non-linear optimization could be used [29][44] and more general lighting models (for 
example including the cone angles from the light sources discussed) could be extracted. It also 
might be worth investigating tying windowing into the fitting process. Initially fit with a 
smoothed version of the function, and then dial back the amount of windowing – effectively 
steering to better local minima. The content sensitive windowing ideas need to be fleshed out, 
particularly when integrating with techniques like PRT. 


One tool that l've found to be invaluable when playing around with Spherical Harmonics is a 
symbolic math program. I’ve used Maple (http://www.maplesoft.com), but other programs like 
Mathematica (http://www.wolfram.com) would work equally well. The DirectX SDK 


(http://msdn.microsoft.com/directx) has functions for evaluation, rotation, products and several 
analytic lighting models, along with samples using both PRT and irradiance environment maps. 
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Appendix A1 Recursive Rules for Evaluating SH Basis Functions 
Recurrence relations [48] can be used to efficiently evaluate the polynomial forms of the SH 
basis functions. Recall the formulas for the basis functions: 


VERE") m>0 (УК cosm RP (cos8) m >o 
yi" = 4 v2Im(¥") т < 0 = AV2k[" біпітір Р," (со80) т < 0 
yo т=0 КОР (cos 6) mes 


Given a point (x,y,z) on the unit sphere, the Associated Legendre polynomials (only depend on Z, 
divided by sin 077) can be evaluated using these recurrences (PE = 1): 
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Where you increment | in the inner loop, m in the outer. 


The trigonometric addition formula can be used to evaluate the dependent terms (multiplied 
by sin Ө" so that they can be computed from (x,y) coordinates on the sphere), where 
5(0) =0,С(0) = 1: 

5(т) = х5(т-1)-усС(т-1) 

С(т) = хС(т-1)-у5(т- 1) 


In the original equation replace sin|m|q with S(m) and cos mq with C(m). It is generally 
more efficient to generate code based on these equations, since it naturally factors common 
subexpressions. 


Appendix А2 Polynomial Forms of SH Basis 
The polynomial forms of the SH basis functions are listed below, L is the band index, M is the 
basis function. Note that Maple randomly changes the order of L and M... 
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Appendix A3 ZH Coefficients for Spherical Light Source 
Given a light source that subtends an angle a in radians here are the symbolic integrals for the 
first 6 bands: 
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Appendix A4 ZH Coefficients for Smooth Cone 

Given a cone that subtends the angle a in radians, the light source has intensity 1 at the north 
pole, and falls off to zero at angle a. At 6™ order this function should not be evaluated with an 
angle less than about 8 degrees using single precision. The derivative of the smoothing function 
is О at the north pole and a. The first 6 bands are: 
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Appendix A5 Solving for Coefficients to Approximate SH 


Environment Map with Directional and Ambient Lights 
You can compute the intensity s that minimizes the approximation error between a directional 
light in direction d and the original lighting environment: 


Е(с) = (Le — cLa)? 


Where L, is the SH representation of the lighting environment, 1,4 is the SH representation of 
the lighting model in direction d. The solution is: 
Le 9 La 
C= 
La ° La 


If you want to add an ambient light, you need to minimize the error function: 


E(e,a) = | (cha(s) ниб) + aLa G) Нуб) — Le(s) * HyG)) ds 


Absorbing the convolution into the lighting and differentiate with respect to each variable: 


13 xcd 1 P " А " 
You can optimize using either radiance (leave vectors alone) or irradiance (convolve vectors.) 
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2-2 | (cha (8) + afa (s) – 1,(5)) L,(s)ds 

Now solve for the two equations equal to zero" to find the minimum. L, (s) is a scalar 


function”, and due to the orthogonality of SH the integrals are all simple dot products of SH 
vectors. This leaves you with the following equations: 


cA+aB =D 
cB+aC=E 
Where: 
5087 16 4 D 205 
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Solving for c and a results in: 
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It turns out there is a simpler way to arrive at the same results. Solve for the intensity of a 
directional light where both the environment and the light don't include the DC term, then 
computing an ambient term that reconstructs the environment DC term when using the scaled 
directional light. This results in: 
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In the above expression the dot product ignores the DC term and Л, [0] is the DC term of the 
environment light. We will use similar techniques for 2 and 3 lights (where the equations get 
much nastier if you handle a at the same time as the light intensities.) 


Multiple lights can be done in a similar manner, first by removing DC from all the light vectors, 
with two lights giving you this error function: 


Есе) = | (colas) + аа) = 1, )) ds 


14 The Hessian matrix is Й 1 which has positive eigenvalues (2 + V2), so this is a minimum. 


* SH vector only has a non-zero in the DC term. 


Then differentiate with respect to each of the coefficients and solve for zero. The following 
expression gives you the intensities: 
A -В 
С] _|A2—B2 А? - В? 22 
C1 —B A Le°Lao 


42-82 Az — Be 
Where: 
А = Lao ° Гао, B = Lao o Lat 


К is worth pointing out that A is constant independent of direction (її depends on order and if 
the vectors are convolved or not.) The ambient term is then: 


a = (L,[0] — CcoLao[0] + алацор) = 


For three lights the intensities are: 
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Where: 
С = Lao ° Laz, В = Lay ° Laz, E = 2BCD + А(А? — B?-C? – D?) 
This matrix is symmetric, and the ambient coefficient is: 


a = (L,[0] — (coLao[0] + c1La1[0] + слао) 


Appendix A6 Least Squares Projection 
First we will show that for orthonormal basis functions, least squares projection actually is 


achieved by simply integrating against the basis functions. We want to find the coefficient 
vector c that minimizes this expression: 


2 


E(c) — ciyi(s) — f(s) | ds 
I» 


1 


This can be done by differentiating with respect to each coefficient, and solving the first 
derivatives for zero? 


dE 
Lee, | 2,0) - f6) |у) ds 


ас, - 
i 


Because the basis functions are orthonormal we know that / y; (s)y; (s)ds = бу exploiting this 
fact and solving for zero we get 


р а | yk Gf (9) ds 


So that direct integration gives you the least squares result. 


Now we derive the coefficient vector g that minimizes an error functional that includes a 
penalty based on the squared Laplacian integrated over the sphere. We know from above the 
coefficient vector c that minimizes pure squared error. We introduce a function h that is used to 
simplify indexing 


[а?а = У Y Pos 120" = Y hf? 
1-0 


l=1 т--1 


The error function then is 


Е(9) = (Se - 23 НЭЭ higi” 


Differentiating you get 


dE 
dui = 2(gxy — ск) + 2Ah, 9k 


Solving for zero you get 


Ck 


9k = x Any) 


Appendix A7 Solving for Lamba to Reduce the Squared Laplacian 


This can be done quite easily using Newtons method. The squared Laplacian is: 


16 . x 5 А 1 z n 1 А Қ 
Second derivative mixed partials are all zero, unmixed partials аге all 1, so the hessian is an identity 
matrix and we have a minima 


1 


n 1 п 
№ = PUHE =) га+1 У ("= LB, 
1=1 1=1 


1=1 m=-1 


3 


т=—1 


The array of L values (Lj = 12(1 + 1)? ) is static, and the array of B values can be computed 


В! = X uh). Newtons method takes an initial guess (0 works well for this problem) and 


refines it using the following recurrence: 


f (a) 


Anti = Àn E Zo) 
7 


Where f is the function we are searching for a root and f’ is the first derivative. Replacing f?” 


with cj" which includes A and factoring you have: 


n 


fA) = 42 — N 
2- (1 + AL)? 
n 2 


3 
£a AL) 


Where A?is the original squared Laplacian. Iterate until a maximum number of iteration has 


occurred, or the absolute value of successive approximations is below a threshold (1e-6 seems 
to work well in practice.) This can be represented as a polynomial (multiplying through by the 
product of all the denominators) of degree (n — 1)?, where n is the order. This would only be 


useful for quadratic and lower orders, where a closed form solution of the roots can be 


computed and even then it is not clear if it would be much faster than an iterative solution. 


Appendix A8 Code for Multiplying SH by Analytic Functions 


Below are two functions, both defined over the hemisphere. The first is simply constant, which 


can be useful in scenes with a dominant ground plane. The second is a clamped cosine in Z. 


// generates 6th order SH coefficients from multiplying 6th order lighting by hemisphere in Z 


void HemiMult(float *R, float *L) 


{ 


R[O] 
R[1] 
R[?] 
R[3] 
R[4] 
R[5] 
R[6] 
R[7] 
R[8] 
R[9] 
R[10] 
R[11] 
R[12] 


Hou pn og go og gp oH dH oH 


0.433012702f * L[2] - 0.1653594569f * 1[12] + 0.1036445247f * L[30] + 0.5f * L[0]; 
9.5f * L[1] + 0.4192627457f * L[5] - 0.1711632992f * L[19]; 

-0.05412658775f * L[20] + 0.5f * L[2] + 0.2420614591f * L[6] + 0.433012702f * L[0]; 
9.5f * L[3] - 0.1711632992f * L[21] + 0.4192627457f * L[7]; 

9.5f * L[4] - 0.1713860232f * L[28] + 0.4133986423f * L[10]; 

-0.06477782793f * L[29] + 0.4192627457f * L[1] + 0.5f * L[5] + 0.2614562582f * L[11]; 
0.58 * L[6] - 0.1448476267f * L[30] + 0.2420614591f * L[2] + 0.3697549864f * L[12]; 
-0.06477782793f * 1[31] + 0.4192627457f * L[3] + 0.5f * L[7] + 0.2614562582f * L[13]; 
9.5f * L[8] - 0.1713860232f * L[32] + 0.4133986423f * L[14]; 

9.5f * L[9] + 0.41015625f * L[17]; 

0.2685102947f * L[18] + 0.5f * L[10] + 0.4133986423f * L[4]; 

0.360244363f * L[19] + 0.5f * L[11] + 0.2614562582f * L[5]; 

0.2790440836f * L[20] + 0.5f * L[12] - 0.1653594569f * L[0] + 0.3697549864f * L[6]; 


} 


R[13] 
R[14] 
R[15] 
R[16] 
R[17] 
R[18] 
R[19] 
R[20] 
R[21] 
R[22] 
R[23] 
R[24] 
R[25] 
R[26] 
R[27] 
R[28] 
R[29] 
R[30] 
R[31] 
R[32] 
R[33] 
R[34] 
R[35] 


ососососоосоосоосооооооосоос 


0.51 
0.51 
0.51 
0.51 


// generates 6th 


.2790440836f * L[12] * 
.2856107251f * L[31] + 
.5f * L[22] + 0.35621916f * L[32] + 0.2685102947f * L[14]; 

.2720668773f * L[33] + 0.5f * L[23] + 0.41015625f * L[15]; 

.5f * L[24] * 0.408100316f * L[34]; 

.5f * L[25]; 

.408100316f * L[16] + 0.5f * L[26]; 

.5f * L[27] + 0.2720668773f * L[17]; 

.5f * L[28] - 0.1713860232f * L[4] + 0.35621916f * L[18]; 

-0.06477782793f * L[5] + 0.5f * L[29] + 0.2856107251f * L[19]; 

0.3498002708f * L[20] + 0.1036445247f * 1[0] - 0.1448476267f * L[6] + 0.5f * L[30]; 
-0.06477782793f * L[7] + 0.2856107251f * L[21] + 0.58 * L[31]; 


.5f * L[13] + 0.360244363f * L[21] + 0.2614562582f * L[7]; 
.2685102947f * L[22] + 0.4133986423f * L[8] + 0.5f * L[14]; 
.41015625f * L[23] * 0.5f * L[15]; 

.408100316f * L[26] + 0.5f * L[16]; 

.2720668773f * L[27] 
.2685102947f * L[10] 
.2856107251f * L[29] - 


+ 0.5f * L[17] + 0.41015625f * L[9]; 

+ 0.35621916f * L[28] + 0.5f * L[18]; 
0.1711632992f * L[1] + 0.5f * L[19] + 0.360244363f * L[11]; 
0.5f * L[20] + 0.3498002708f * L[30] - 0.05412658775f * L[2]; 
0.360244363f * L[13] + 0.5f * L[21] - 0.1711632992f * L[3]; 


* L[32] - 0.1713860232f * L[8] + 0.35621916f * L[22]; 
* L[33] + 0.2720668773f * L[23]; 

* L[34] + 0.408100316f * L[24]; 

* L[35]; 


order SH coefficients from multiplying 6th order lighting by clamped cos in Z 


void CosMult(float *R, float *L) 


{ 


// constants to make code fit 


const 
const 
const 
const 
const 
const 
const 
const 
const 
const 
const 
const 
const 


R[O] 
R[1] 
R[2] 
R[3] 
R[4] 
R[5] 
R[6] 
R[7] 
R[8] 
R[9] 
R[10] 
R[11] 
R[12] 
R[13] 
R[14] 
R[15] 
R[16] 
R[17] 
R[18] 
R[19] 
R[29] 
R[21] 
R[22] 
R[23] 
R[24] 
R[25] 
R[26] 
R[27] 


float T1 = 0.09547032698f; 

float T2 - 0.1169267933f; 

float T3 - 0.2581988897f; 

float T4 = 0.2886751347f; 

float T5 = 0.2390457218f; 

float T6 = 0.2535462764f; 

float T7 = 0.2182178903f; 

float T8 - 0.1083940384f; 

float T9 = 0.2519763153f; 

float TA = 0.2439750183f; 

float TB = 0.3115234375f; 

float TC = 0.2512594538f; 

float TD = 0.31640625f; 

0.25f * L[0] - 0.03125f * L[20] + ТА * L[2] + 0.1397542486f * L[6]; 
0.2236067977f * L[5] + T2* L[11] - 0.02896952533f * L[29] + 0.1875f * L[1]; 
90.375f * L[2] + T4 * L[0] + T3 * L[6] - 0.01495979856f * L[30] + T1 * L[12]; 
0.2236067977f * L[7] - 0.02896952533f * L[31] + T2* L[13] + 0.1875f * L[3]; 
0.101487352f * L[18] + 0.1889822365f * L[10] + 0.15625f * L[4]; 

9.2236067977Ғ * L[1] + T5 * L[11] + 0.3125f * L[5] + 0.09568319309f * L[19]; 
T6 * L[12] + T3 * L[2] + 8.3125f * L[6] + 0.113550327f * L[20] + 0.1397542486f * L[0]; 
T5 * L[13] + 0.09568319309f * L[21] + 0.2236067977f * L[3] + 0.3125f * L[7]; 
9.15625Ғ * L[8] + 0.101487352f * L[22] + 0.1889822365f * L[14]; 

0.09068895910f * L[27] + 0.1666666667f * L[17] + 0.13671875f * L[9]; 

= T7 * L[18] + 0.2734375f * L[10] + 0.09068895910f * L[28] + 0.1889822365f * L[4]; 
= T8 * L[29] + TA * L[19] + T5 * L[5] + 0.30078125f * L[11] + T2* L[1]; 

= T9 * L[20] + 0.328125f * L[12] + T6 * L[6] + 0.1028316139f * L[30] + T1 * L[2]; 

= 0.30078125f * L[13] + T8 * L[31] + TA * L[21] + T5 * L[7] + T2* L[3]; 

= 0.2734375f * L[14] + 0.09068895910f * L[32] + 0.1889822365f * L[8] + T7 * L[22]; 

= 0.0906889591Ғ * L[33] + 0.13671875f * L[15] + 0.1666666667f * L[23]; 

= 0.1507556723f * L[26] + 0.123046875f * L[16]; 

- 0.24609375f * L[17] * 0.1666666667f * L[9] * 0.201007563f * L[27]; 

= 0.101487352f * L[4] + T7 * L[10] + 0.2302830933f * L[28] + 0.28125f * L[18]; 
= TD * L[19] + 0.09568319309f * L[5] + TA * L[11] + 0.2461829819f * L[29]; 

= 79 * L[12] + TC * L[30] + 0.113550327f * L[6] - 0.03125f * L[O] + TD * L[20]; 
- 0.2461829819f * L[31] * TA * L[13] * 0.9568319309e-1f * L[7] * TD * L[21]; 

= 0.101487352f * L[8] + 0.28125f * L[22] + T7 * L[14] + 0.2302830933f * L[32]; 
= 0.201007563f * L[33] + 0.1666666667f * L[15] + 0.24609375f * L[23]; 

= 0.1507556723f * L[34] + 0.123046875f * L[24]; 

- 0.1127929688f * L[25]; 

- 0.2255859375f * L[26] * 0.1507556723f * L[16]; 

= 0.9068895910е-1Ғ * L[9] + 0.201007563f * L[17] + 0.2631835938f * L[27]; 


R[28] 
R[29] 
R[30] 
R[31] 
R[32] 
R[33] 
R[34] 
R[35] 


HoH Hn Ho dg gb gp Ogg 


оооооооо 


.2302830933f * L[18] + 0.9068895910e-1f * L[10] + 0.30078125f * L[28]; 
.1083940385f * 1[11] + 0.2461829819f * L[19] - 0.02896952533f * L[1] + ТВ * L[29]; 
.322265625Ғ * L[30] + TC * L[20] - 0.01495979856f * L[2] + 0.1028316139f * L[12]; 
.2461829819f * L[21] - 0.02896952533f * L[3] + TB * L[31] + 0.1083940385f * L[13]; 
.09068895910f * L[14] + 0.2302830933f * L[22] + 0.30078125f * L[32]; 

.09068895910f * L[15] + 0.201007563f * L[23] + 0.2631835938f * L[33]; 

.2255859375f * L[34] * 0.1507556723f * L[24]; 

.1127929688f * L[35]; 


Appendix A9 Ambient Cube Basis 


The Ambient Cube basis is used by Valve [26]; it consists of 6 basis functions, each defined over 


a hemisphere: 


2 
" x»0 V1= 


2. E x0 


x? x<0 


2 82550 
v2= {у У>0,уз= | У 
0 у<0 ус у<0 


2 
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This basis is not orthogonal, so coefficients cannot be generated by simply integrating against 


the basis functions. To compute optimal projection coefficients you solve a linear least squares 


problem (similar to Appendix A6 Least Squares Projection, but without orthonormal basis 


functions): 


2 


во = | [au - r6 | as 
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Where the V; (s) are the Ambient Cube basis functions. Differentiating you get: 


сЕ | У ails) - fG) И, (5) ds 


Integration is linear so you can re-arange the terms and solve for zero getting 


Уа иом) = | Gv. Gods 


The left hand side gives you a row of a matrix A, where Aj; = BAQ (s)ds and the right hand 
side will give you a vector (integrals of function times this basis function.) This leads to a linear 
system: Ac = b the inverse of A is: 


11 1 3 3 3 3 
4n 4n 8л 8л 8л 8л 
1 11 З З З З 
4n 4n — Om 8n 8л 8л 
З З 11 1 З З 
8л 8л 4m 4т 8n 8n 
3 3 1 11 3 3 
КҮ: 8л 4m 4n 8л 8m 
3 3 3 3 11 1 
8n 8n 8л 8л 4т 4n 
3 3 3 3 1 11 
8n 8n 8л 8n 4т 4n 


To project a function represented using this basis into SH you just integrate the V's against the 
SH basis, for quadratic spherical harmonics that results in the following matrix: 


d 


JE 


3 3 ЕЕ 3 3 E 
З үл 3 4T 
0 0 2 А 0 0 
0 0 0 0 (3 4n (3 үт 
4 4 
Зуп 3 4T 
à л 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
л 45 л 45 dm 4/5 іл(Б  24m4[5 п 45 
15 15 15 15 15 15 
0 0 0 0 0 0 
(іт 4/15 л 4/15 п/15 4m415 1 7 
15 15 15 15 


The even degrees above quadratic are т the null-space of the Ambient Cube basis, and the odd 
degrees above quadratic are in the null space of the clamped cosine function, so if irradiance 
environment maps are being used, no higher order is necessary. The Ambient Cube basis can 
exactly reconstruct the DC term, approximate the linear term, exactly reconstructs two of the 
quadratic basis functions (y2, 92) and has three of the quadratic basis functions in the null space 


(752,21, yd). 


To project from SH into the Ambient Cube basis you need to minimize this error function: 


2 


во = | [Savi - Уу | as 
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Where | are the SH lighting coefficients (assumed to be convolved.) Once again differentiate 


with respect to the unkown variables: 


dn 2 | 2,440) “2H ү, (5) ds 


i 


And then solve for zero: 


Уа [uec = У | vos 
i j 


This results in a linear system as well, Ac = ВІ, where A is the same as before B is a matrix 
where By = Ли. ($); (s)ds. The matrix 4-1 B can be used to move from SH to the Ambient 


Cube basis: 
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Appendix A10 Shader/CPU code for Irradiance Environment Maps 
Given a quadratic SH representation of a lighting environment, it is fairly simple to generate 
shader code. In [35] a matrix representation is used, however this turns out to require more 
instructions" (15 vs. 11) and more constants (12 vs. 7) compared to a more direct evaluation of 
the quadratic spherical harmonics. You fold leading constants of the polynomials into the 
lighting coefficients, group all but ys" by channel (using float4's) and keep ys* as a color and 


On scalar GPU's this gap is even more significant — 60 vs. 42 


fold part of у? into DC. Shader code for evaluation and CPU code to setup the constants is 
shown below. The normal passed in to shader code should be normalized and the 4^ channel 
should be 1.0. The CPU code takes an array of 3 float pointers to quadratic radiance SH 
coefficients and an Effect to bind the constants to. 


// constants containing irradiance environmne map 
float4 cAr; 
float4 cAg; 
float4 cAb; 
float4 cBr; 
float4 cBg; 
float4 cBb; 
float4 cC; 


float3 ShadeIrad(float4 vNormal) 


{ 
float3 x1, x2, x3; 


// Linear + constant polynomial terms 
х1.г = dot(cAr,vNormal); 
X1.g = dot(cAg, vNormal); 
X1.b = dot(cAb, vNormal); 


// 4 of the quadratic polynomials 
float4 vB - vNormal.xyzz * vNormal.yzzx; 
x2.r = dot(cBr,vB); 

X2.g = dot(cBg, vB); 

х2.0 = dot(cBb, vB); 


// Final quadratic polynomial 

float vC - vNormal.x*vNormal.x - vNormal.y*vNormal.y; 
X3 - cC.rgb * vC; 

return х1+х2+х3; 


void SetSHEMapConstants( float* fLight[3], ID3DXEffect* pEffect) 
{ 

// Lighting environment coefficients 

D3DXVECTORA vCoeff[3]; 


static const float s fSqrtPI - ((float)sqrtf(D3DX PI)); 
const float fCO 1.0f/(2.0f*s fSqrtPI); 

const float fC1 (float)sqrt(3.0f)/(3.0f*s fSqrtPI); 
const float fC2 (float)sqrt(15.0f)/(8.0f*s fSqrtPI); 
const float fC3 (float)sqrt(5.0f)/(16.0f*s fSqrtPI); 


const float ТС4 9.5f*fC2; 
int iC; 
for( 1С-0; iC«3; 1С44 ) 
{ 
vCoeff[iC].x = -fC1*fLight[ic][3]; 
vCoeff[ic].y = -fC1*fLight[ic][1]; 
vCoeff[iC].z = fC1*fLight[ic][2]; 
vCoeff[ic].w = fco*fright[ic][0] - fC3*fLight[ic][6]; 


} 

pEffect->SetVector( "cAr", &vCoeff[0] ); 
pEffect->SetVector( "cAg", &vCoeff[1] ); 
pEffect->SetVector( "cAb", &vCoeff[2] ); 


for( 1С-0; iC<3; iC++ ) 


{ 
vCoeff[iC].x = fc2*fLight[iC][4]; 
vCoeff[iC].y - -fc2*fLight[ic][5]; 
vCoeff[iC].z = 3.0f*fC3*fLight[iC][6]; 
vCoeff[iC].w - -fc2*fLight[ic][7]; 
j 


pEffect->SetVector( "cBr", &vCoeff[0] ); 
pEffect->SetVector( "cBg", &vCoeff[1] ); 
pEffect->SetVector( "cBb", &vCoeff[2] ); 


vCoeff[0].x = fc4*fLight[0][8]; 
vCoeff[0].y = fCA*fLight[1][8]; 
vCoeff[0].z = fC4*fLight[2][8]; 
vCoeff[0].w = 1.0f; 


pEffect->SetVector( "cC", &vCoeff[0] ); 


Version History 


e 2/14/2008 first public version 

e 4/1/2008 Added appendix with recurrence relations when evaluating SH, renumbered appendices 
e 3/17/2009 Fixed typos in the recurrence relations 

e 11/30/2009 Ambient light and optimal linear direction typos fixed 


e 2/10/2010 Typo in defining product matrix (thanks to Don Holden for pointing it out) 


